function [vadist,ll,gr,scr] = emstep(vadist,Y,H,Q)

[numvars,numcomp] = size(vadist.mean);
n = numel(Y);

E = zeros([numvars n numcomp]);
V = zeros([numvars numvars n numcomp]);
like = zeros([n numcomp]);
for c = 1:numcomp
    MU = vadist.mean(:,c);
    SIGMA = vadist.cov(:,:,c);
    for i = 1:n
        [E(:,i,c),V(:,:,i,c),like(i,c)] = mvnfilter(Y{i},H{i},MU,SIGMA,Q{i});
    end
end
like = bsxfun(@times,like,vadist.pr);
qi = bsxfun(@rdivide,like,sum(like,2));
ll = log(sum(like,2));

ll = sum(ll);
E2 = V + bsxfun(@times,reshape(E,[numvars 1 n numcomp]),reshape(E,[1 numvars n numcomp]));

if nargout>=3

    gr_mean = zeros([numvars numcomp n]);
    gr_cov = zeros([numvars numvars numcomp n]);
    for c = 1:numcomp
        MU = vadist.mean(:,c);
        SIGMA = vadist.cov(:,:,c);
        SIGMAinv = SIGMA \ eye(size(SIGMA));
        gr_mean(:,c,:) = reshape(SIGMAinv*bsxfun(@minus,E(:,:,c),MU),[numvars 1 n]);
        for i = 1:n
            eC = E2(:,:,i,c) - E(:,i,c)*MU' - MU*E(:,i,c)' + MU*MU';
            gr_cov(:,:,c,i) = -.5*(SIGMAinv - SIGMAinv*eC*SIGMAinv);
        end
    end
    gr_mean = bsxfun(@times,reshape(qi',[1 numcomp n]),gr_mean);
    gr_cov = bsxfun(@times,reshape(qi',[1 1 numcomp n]),gr_cov);

    if nargout==4

        scr_mean = reshape(gr_mean,[numvars*numcomp n]);
        scr_cov = reshape(gr_cov,[numvars^2 numcomp n]);
        scr_cov(triu(true(numvars),1),:,:) = 2*scr_cov(triu(true(numvars),1),:,:);
        scr_cov = reshape(scr_cov(triu(true(numvars)),:,:),[numcomp*(numvars^2 + numvars)/2 n]);

        scr_pr = bsxfun(@rdivide,qi,vadist.pr);
        scr_pr = bsxfun(@minus,scr_pr(:,1:end-1),scr_pr(:,end));

        scr = [scr_mean' scr_cov' scr_pr];

    end

    gr_mean = sum(gr_mean,3);
    gr_cov = sum(gr_cov,4);
    gr_cov_chol = cell([1 numcomp]);
    for c = 1:numcomp
        SIGMA_chol = chol(vadist.cov(:,:,c));
        gr_cov_chol{c} = 2*SIGMA_chol*gr_cov(:,:,c);
        gr_cov_chol{c}(eye(numvars)==1) = diag(gr_cov_chol{c}).*diag(SIGMA_chol);
    end
    gr_cov_chol = cat(3,gr_cov_chol{:});
    
    gr_pr_log = sum(bsxfun(@minus,qi(:,1:end-1),vadist.pr(1:end-1)),1);

    gr = [gr_mean(:); 
          gr_cov_chol(repmat(triu(true(numvars)),[1 1 numcomp]));
          gr_pr_log'];
    
end

vadist.pr = mean(qi,1);
for k = 1:numcomp
    vadist.mean(:,k) = (E(:,:,k)*qi(:,k))./sum(qi(:,k));
end
for k = 1:numcomp
    w_k = qi(:,k);
    E2_k = reshape(reshape(E2(:,:,:,k),[numvars^2 n])*w_k,[numvars numvars])./sum(w_k);
    E1E0_k = ((E(:,:,k)*w_k)./sum(w_k))*vadist.mean(:,k)';
    E0E0_k = vadist.mean(:,k)*vadist.mean(:,k)';
    vadist.cov(:,:,k) = E2_k - E1E0_k - E1E0_k' + E0E0_k;
    vadist.cov(:,:,k) = (vadist.cov(:,:,k) + vadist.cov(:,:,k)')./2;
end
